
module grmult
!-----------------------------------------------------------------------
!
! Purpose:
! Compute non-linear dynamics terms in grid point space (in preparation
! for SLD interpolation). This is broken up into two interfaces an
! initialization phase for constants that will be used and the run phase
! where the non-linear terms are actually computed.
!
!---------------------------Code history--------------------------------
!
! Author:  J. Olson
!-----------------------------------------------------------------------
!
! $Id$
! $Author$
!
!-----------------------------------------------------------------------

   use shr_kind_mod,    only: r8 => shr_kind_r8
   use pmgrid,          only: plon, plev, plevp
   use physconst,       only: rair
   use hycoef,          only : hyam, hybm, ps0, hyai, hybi, hypd, hybd, hypm, hypi, nprlev
   implicit none

   private
!
! Public interfaces
!
   public grmult_init   ! Initialization
   public grmult_run    ! Run method
!
! Private module data
!
   real(r8) :: hortalc(plev)           ! analytic Hortal temperature correction term
   real(r8) :: hdel   (0:plev-1)       ! del-"hortalc" in vertical
   real(r8) :: onemeps                 ! 1 - epssld (SLD decentering coefficient)
   real(r8) :: onepeps                 ! 1 + epssld (SLD decentering coefficient)
   real(r8) :: detai(plevp)            ! interval between interfaces
   real(r8) :: facm1                   ! interpolation factor for time n-1
   real(r8) :: facm2                   ! interpolation factor for time n-2

contains

subroutine grmult_init()
!-----------------------------------------------------------------------
!
! Purpose:
! Initialization of constants used for grmult.
!
!-----------------------------------------------------------------------
   use physconst, only: rga
   use scanslt,   only: epssld

!-----------------------------------------------------------------------
!
!  Local variables
!
!
   real(r8) tmp            ! tmp variable
   real(r8) tmp1           ! tmp variable
   real(r8) psref          ! reference surf pres in Hortal calc.
   real(r8) pref           ! reference pres at level k in Hortal calc
   real(r8) tnot           ! reference surf temp in Hortal calc.
   real(r8) tref           ! reference temp at level k in Hortal calc.
   real(r8) lapse          ! standard atmosphere lapse rate
   real(r8) tcut           ! cutoff temperature below which "Hortal"
   integer :: k            ! Index
!
! Initialize array of constants for Hortal temperature correction
!
   psref = 101325._r8
   tnot  = 288._r8
   lapse = 0.0065_r8
   tcut  = 216.5_r8
   tmp   = rair*lapse*rga
!
   do k = 1,plev
      pref = hyam(k)*ps0 + hybm(k)*psref
      tmp1 = pref/psref
      tref = tnot*tmp1**tmp
      if (tref > tcut) then
         hortalc(k) = hybm(k)*tmp*tref/tmp1
      else
         hortalc(k) = 0._r8
      endif
   end do
!
! Compute vertical difference
!
   hdel(0) = hortalc(1)
   do k = 1,plev-1
      hdel(k) = hortalc(k+1) - hortalc(k)
   end do
   onemeps = 1._r8 - epssld
   onepeps = 1._r8 + epssld
   facm1   =  3._r8/2._r8
   facm2   = -1._r8/2._r8
!
   do k = 1,plev
      detai  (k) = (hyai(k+1) + hybi(k+1)) - (hyai(k) + hybi(k))
   end do

end subroutine grmult_init

!
!-----------------------------------------------------------------------
!

subroutine grmult_run (ztodt   ,rcoslat ,coslat  ,coriol  ,div     , &
                   q3      ,t3      ,u3      ,v3      , &
                   ed1     ,ql      ,qm      ,tl      , &
                   tm      ,phis    ,phisl   ,phism   ,dpsl    , &
                   dpsm    ,omga    ,pmid    ,pdel    ,ps      , &
                   logps   ,rpmid   ,etamid  ,fu      ,fv      , &
                   t2      ,t3m1    ,u3m1    ,v3m1    ,etadot  , &
                   dpslon  ,dpslat  ,urhs    ,vrhs    ,trhs    , &
                   prhs    ,lnpssld ,prhssld ,tarrsld , &
                   parrsld ,ktoop   ,nlon    )
!-----------------------------------------------------------------------
!
! Purpose:
! Compute non-linear dynamics terms in grid point space (in preparation
! for SLD interpolation)
!
!
!-----------------------------------------------------------------------
  use physconst,    only: zvir, cappa, cpvir
  use time_manager, only: is_first_step
  use commap,       only: bps, t0, href


!------------------------------Arguments--------------------------------
!
  real(r8), intent(in)   :: ztodt                 ! delta-t
  real(r8), intent(in)   :: rcoslat               ! 1/(cos lat)
  real(r8), intent(in)   :: coslat                ! cos(lat)
  real(r8), intent(in)   :: coriol                ! Coriolis parameter
  real(r8), intent(in)   :: div     (plon,plev)   ! divergence
  real(r8), intent(in)   :: q3      (plon,plev)   ! specific humidity
  real(r8), intent(in)   :: t3      (plon,plev)   ! temperature
  real(r8), intent(in)   :: u3      (plon,plev)   ! zonal wind
  real(r8), intent(in)   :: v3      (plon,plev)   ! meridional wind
  real(r8), intent(in)   :: ed1     (plon,plev)   ! (1/ps)etadot(dp/deta) (time n-1)
  real(r8), intent(in)   :: ql      (plon,plev)   ! longitudinal derivative of q
  real(r8), intent(in)   :: qm      (plon,plev)   ! latitudinal  derivative of q
  real(r8), intent(in)   :: tl      (plon,plev)   ! zonal derivative of T
  real(r8), intent(in)   :: tm      (plon,plev)   ! meridional derivative of T
  real(r8), intent(in)   :: phis    (plon)        ! Phi at surface
  real(r8), intent(in)   :: phisl   (plon)        ! longitudinal derivative of phis
  real(r8), intent(in)   :: phism   (plon)        ! latitudinal  devivative of phis
  real(r8), intent(in)   :: dpsl    (plon)        ! longitudinal component of grad ln(ps)
  real(r8), intent(in)   :: dpsm    (plon)        ! latitudinal  component of grad ln(ps)
  real(r8), intent(in)   :: omga    (plon,plev)   ! vertical pressure velocity
  real(r8), intent(in)   :: pmid    (plon,plev)   ! pressure at full levels
  real(r8), intent(in)   :: pdel    (plon,plev)   ! layer thicknesses (pressure)
  real(r8), intent(in)   :: ps      (plon)        ! Surface pressure (n)
  real(r8), intent(in)   :: logps   (plon)        ! log(ps)
  real(r8), intent(in)   :: rpmid   (plon,plev)   ! 1./pmid
  real(r8), intent(in)   :: etamid  (plev)        ! midpoint values of eta (a+b)
  real(r8), intent(inout):: fu      (plon,plev)   ! nonlinear term - u momentum eqn
  real(r8), intent(inout):: fv      (plon,plev)   ! nonlinear term - v momentum eqn
  real(r8), intent(inout):: t2      (plon,plev)   ! nonlinear term - temperature
  real(r8), intent(out)  :: t3m1    (plon,plev)   ! T at time n-1
  real(r8), intent(inout):: u3m1    (plon,plev)   ! U at previous time step
  real(r8), intent(inout):: v3m1    (plon,plev)   ! V at previous time step
  real(r8), intent(out)  :: etadot  (plon,plevp)  ! vertical velocity in eta coordinates
  real(r8), intent(out)  :: dpslon  (plon,plev)   ! Pressure gradient term
  real(r8), intent(out)  :: dpslat  (plon,plev)   ! Pressure gradient term
  real(r8), intent(inout):: urhs    (plon,plev)   ! RHS of U  eqn valid for mid-point
  real(r8), intent(inout):: vrhs    (plon,plev)   ! RHS of V  eqn valid for mid-point
  real(r8), intent(inout):: trhs    (plon,plev)   ! RHS of T  eqn valid for mid-point
  real(r8), intent(inout):: prhs    (plon,plev)   ! RHS of Ps eqn valid for mid-point
  real(r8), intent(out)  :: lnpssld (plon,plev)   ! RHS Ps term for SLD
  real(r8), intent(out)  :: prhssld (plon,plev)   ! RHS Ps term for SLD
  real(r8), intent(out)  :: tarrsld (plon,plev)   ! T  at arr. pt. (SLD)
  real(r8), intent(out)  :: parrsld (plon,plev)   ! Ps at arr. pt. (SLD)
  real(r8), intent(out)  :: ktoop   (plon,plev)   ! (Kappa*T)*(omega/P)
  integer , intent(in)   :: nlon                  ! number of longitudes for this latitude
!
!---------------------------Local workspace-----------------------------
!
  real(r8) tmp1                  ! temporary workspace
  real(r8) tmp2                  ! temporary workspace
  real(r8) tmp                   ! temporary workspace
  real(r8) tmpk                  ! workspace
  real(r8) tmpkp1                ! workspace
  real(r8) u3l     (plon,plev)   ! u-wind used locally only
  real(r8) v3l     (plon,plev)   ! v-wind used locally only
  real(r8) tv      (plon,plev)   ! virtual temperature
  real(r8) ddpk    (plon)        ! partial sum of div*delta p
  real(r8) ddpn    (plon)        ! complete sum of div*delta p
  real(r8) vkdp    (plon,plev)   ! V dot grad(ln(ps))
  real(r8) vpdsk   (plon)        ! partial sum  V dot grad(ln(ps)) delta b
  real(r8) vpdsn   (plon)        ! complete sum V dot grad(ln(ps)) delta b
  real(r8) lpsstar (plon)        ! Reference ln(Ps) (used to define a new 
!                                ! perturbation Ps)
  real(r8) lpsstarl(plon)        ! long. grad of reference ln(Ps)
  real(r8) lpsstarm(plon)        ! lat.  grad of reference ln(Ps)
  real(r8) rtv     (plon,plev)   ! rair*(tv+t0)
  real(r8) dt                    ! time step
  real(r8) hsl     (plon,plev)   ! zonal      deriv of hydrostatic term
  real(r8) hsm     (plon,plev)   ! meridional deriv of hydrostatic term
  real(r8) pspsl   (plon)        ! Ps*d(lnPs)/d(long.)
  real(r8) pspsm   (plon)        ! Ps*d(lnPs)/d(lat. )
  real(r8) ed1p    (plon,plevp)  ! (1/ps)etadot(dp/deta) (time n-1)
  real(r8) rtvl                  ! zonal      derivative of R*Tv
  real(r8) rtvm                  ! meridional derivative of R*Tv
  real(r8) abp0                  ! constant for grad(H(n)) matrix
!
! Arrays which hold results from the RHS of the prognostic equations.
! Results will be interpolated to trajectory departure points in the routine
! SCANSLT.
!
  real(r8) tsld0a  (plon,plev)   ! RHS of T  eqn valid for mid-point
  real(r8) usldm   (plon,plev)   ! RHS of U  eqn valid for departure pt (n)
  real(r8) vsldm   (plon,plev)   ! RHS of V  eqn valid for departure pt (n)
  real(r8) tsldm   (plon,plev)   ! RHS of T  eqn valid for departure pt (n)
  real(r8) psldm   (plon,plev)   ! RHS of Ps eqn valid for departure pt (n)
  real(r8) urhsl   (plon,plev)   ! RHS of U  eqn valid for mid-point (n+1/2)
  real(r8) vrhsl   (plon,plev)   ! RHS of V  eqn valid for mid-point (n+1/2)
  real(r8) trhsl   (plon,plev)   ! RHS of T  eqn valid for mid-point (n+1/2)
  real(r8) prhsl   (plon,plev)   ! RHS of Ps eqn valid for mid-point (n+1/2)


  integer i,l,k                  ! longitude, level indices
  integer npr                    ! index
!
!-----------------------------------------------------------------------
!
!
! Compute U/V for time n + 1/2
! This will be used locally in the Hortal Temperature correction
!
  do k = 1,plev
     do i = 1,nlon
        u3l  (i,k) = facm1*u3(i,k) + facm2*u3m1(i,k)
        v3l  (i,k) = facm1*v3(i,k) + facm2*v3m1(i,k)
     end do
  end do
!
! Zero auxiliary fields
!
  tmp = 1._r8/(rair*t0(plev))
  do i=1,nlon
     ddpk    (i)     = 0._r8
     ddpn    (i)     = 0._r8
     vpdsk   (i)     = 0._r8
     vpdsn   (i)     = 0._r8
     lpsstar (i)     = -phis (i)*tmp
     lpsstarl(i)     = -phisl(i)*tmp
     lpsstarm(i)     = -phism(i)*tmp
     pspsl   (i)     = ps(i)*dpsl(i)
     pspsm   (i)     = ps(i)*dpsm(i)
     etadot(i,1)     = 0._r8
     ed1p  (i,1)     = 0._r8
     etadot(i,plevp) = 0._r8
  end do
!
! Virtual temperature
!
  call virtem(nlon, plon, plev, t3      ,q3      ,zvir    ,tv)
!
! calculate some auxiliary quantities
!
  do k=1,plev
     do i=1,nlon
        ed1p(i,k+1) = ed1(i,k)
        rtv(i,k) = rair*tv(i,k)
!
! sum(plev)(div(k)*dp(k))
!
        ddpn(i) = ddpn(i) + div(i,k)*pdel(i,k)
     end do
  end do
!
! sum(plev)(v(k)*grad(lnps)*db(k))
!
  do k=nprlev,plev
     do i=1,nlon
        vkdp(i,k)= rcoslat*(u3(i,k)*pspsl(i) + v3(i,k)*pspsm(i))
        vpdsn(i) = vpdsn(i) + vkdp(i,k)*hybd(k)
     end do
  end do
!
! Compute etadot (top and bottom = 0.)
!
  do k = 1,plev-1
!
! Compute etadot(dp/deta)(k+1/2) and sum(k)(div(j)*dp(j))
!
     do i=1,nlon
        ddpk(i) = ddpk(i) + div(i,k)*pdel(i,k)
        etadot(i,k+1) = -ddpk(i)
     end do
!
! sum(k)(v(j)*grad(ps)*db(j))
!
     if (k.ge.nprlev) then
        do i=1,nlon
           vpdsk(i) = vpdsk(i) + vkdp(i,k)*hybd(k)
           etadot(i,k+1) = etadot(i,k+1) - vpdsk(i) + hybi(k+1)*(ddpn(i)+vpdsn(i))
        end do
     end if
!
! Convert eta-dot(dp/deta) to eta-dot
!
     tmp = etamid(k+1) - etamid(k)
     do i = 1,nlon
        etadot(i,k+1) = etadot(i,k+1)*tmp/(pmid(i,k+1) - pmid(i,k))
     end do
  end do
!
! Zonal and meridional derivatives of the hydrostatic term in the
! momentum equations.
!
  do k = 1,plev
     do i = 1,nlon
        tmp1     = (1._r8 + zvir*q3(i,k))
        tmp2     = t3(i,k)*zvir
        rtvl     = rair*( tmp1*tl(i,k) + tmp2*ql(i,k) )
        rtvm     = rair*( tmp1*tm(i,k) + tmp2*qm(i,k) )
!
        tmp      = rpmid(i,k)*pdel(i,k)
        hsl(i,k) = -rtvl*tmp
        hsm(i,k) = -rtvm*tmp
     end do
     if(k .ge. nprlev) then
        abp0 = ps0*(hyam(k)*(hybi(k+1) - hybi(k)) - hybm(k)*(hyai(k+1) - hyai(k)))
        do i = 1,nlon
           tmp      = rtv(i,k)*abp0/(pmid(i,k)*pmid(i,k))
           hsl(i,k) = hsl(i,k) - pspsl(i)*tmp
           hsm(i,k) = hsm(i,k) - pspsm(i)*tmp
        end do
     endif
  end do
!
! Calculate RHS of all prognostic eqns
!
  do k = 1,plev
     tmp  = cappa*t0(k)*hypi(plevp)/hypm(k)
     tmp1 = hypd(k)/hypi(plevp)
     do i = 1,nlon
!
! Surface Pressure eqn
!
        prhsl(i,k) =  div(i,k)*(tmp1 - pdel(i,k)/ps(i))
        psldm(i,k) = -div(i,k)*tmp1 - ed1p(i,k+1) + ed1p(i,k)
!
! Temperature eqn
!
        ktoop (i,k) = cappa*tv(i,k)/(1._r8 + cpvir*q3(i,k))*omga(i,k)*rpmid(i,k)
        trhsl (i,k) = ktoop(i,k)
        trhsl (i,k) = trhsl(i,k) - omga(i,k)/ps(i)*tmp
        tsldm (i,k) = 0.5_r8*(ed1p(i,k+1) + ed1p(i,k))*tmp
!
! ... horizontal advection portion of Hortal Temperature correction
!
        trhsl (i,k) = trhsl(i,k) - &
                      rcoslat*( u3(i,k)*lpsstarl(i) + v3(i,k)*lpsstarm(i) )*hortalc(k)
!
! ... Ritchie damping term for Temperature eqn.
!
        tsld0a(i,k) = rcoslat*( u3l(i,k)*lpsstarl(i) + v3l(i,k)*lpsstarm(i) )
!
! U/V eqns (includes only the diagonal portion of the hydrostatic term)
!
        urhsl(i,k) = 0.5_r8*hsl(i,k) + href(k,k)*tl(i,k) + bps(k)*dpsl(i)
        vrhsl(i,k) = 0.5_r8*hsm(i,k) + href(k,k)*tm(i,k) + bps(k)*dpsm(i)
        usldm(i,k) = -phisl(i) + v3(i,k)*coriol*coslat -href(k,k)*tl(i,k) - bps(k)*dpsl(i)
        vsldm(i,k) = -phism(i) - u3(i,k)*coriol*coslat -href(k,k)*tm(i,k) - bps(k)*dpsm(i)
     end do
!
! Add pressure gradient terms to momentum tendencies
!
     if (k.ge.nprlev) then
        do i=1,nlon
           tmp        = rtv(i,k)*hybm(k)*rpmid(i,k)
           dpslon(i,k) = rcoslat*tmp*pspsl(i)
           dpslat(i,k) = rcoslat*tmp*pspsm(i)
           urhsl(i,k) = urhsl(i,k) - dpslon(i,k)*coslat
           vrhsl(i,k) = vrhsl(i,k) - dpslat(i,k)*coslat
        end do
     else
        do i = 1,nlon
           dpslon(i,k) = 0._r8
           dpslat(i,k) = 0._r8
        end do
     end if
  end do
!
! Interior levels of the hydrostatic term
!
  do k=1,plev-1
     do l=k+1,plev
        do i=1,nlon
           urhsl(i,k) = urhsl(i,k) + href(l,k)*tl(i,l) + hsl(i,l)
           vrhsl(i,k) = vrhsl(i,k) + href(l,k)*tm(i,l) + hsm(i,l)
!
           usldm(i,k) = usldm(i,k) - href(l,k)*tl(i,l)
           vsldm(i,k) = vsldm(i,k) - href(l,k)*tm(i,l)
        end do
     end do
  end do
!
! Compute vertical advection portion of Hortal Temperature correction
!
  npr = nprlev - 1
  if(npr .lt. 1) npr = 1
  do k = npr,plev-1
     tmpk   = 0.5_r8*hdel(k-1)/detai(k)
     tmpkp1 = 0.5_r8*hdel(k  )/detai(k)
     do i = 1,nlon
        trhsl(i,k) = trhsl(i,k) - (etadot(i,k  )*tmpk + etadot(i,k+1)*tmpkp1)*lpsstar(i)
     end do
  end do
!
! ... bottom level
!
  tmpk = 0.5_r8*hdel(plev-1)/detai(plev)
  do i = 1,nlon
     trhsl(i,plev) = trhsl(i,plev) - etadot(i,plev)*tmpk*lpsstar(i)
  end do
!
! Compute (by extrapolation) RHS terms for time n + 1/2
!
  if(is_first_step()) then
     do k = 1,plev
        do i = 1,nlon
           urhs  (i,k) = urhsl (i,k)
           vrhs  (i,k) = vrhsl (i,k)
           trhs  (i,k) = trhsl (i,k)
           prhs  (i,k) = prhsl (i,k)
        end do
     end do
  endif
!
  do k = 1,plev
     do i = 1,nlon
        tmp         =       urhsl (i,k)
        urhsl (i,k) = facm1*urhsl (i,k) + facm2*urhs  (i,k)
        urhs  (i,k) = tmp
        tmp         =       vrhsl (i,k)
        vrhsl (i,k) = facm1*vrhsl (i,k) + facm2*vrhs  (i,k)
        vrhs  (i,k) = tmp
        tmp         =       trhsl (i,k)
        trhsl (i,k) = facm1*trhsl (i,k) + facm2*trhs  (i,k)
        trhs  (i,k) = tmp
        tmp         =       prhsl (i,k)
        prhsl (i,k) = facm1*prhsl (i,k) + facm2*prhs  (i,k)
        prhs  (i,k) = tmp
     end do
  end do
!
! Prepare SLD arrays for interpolation by SCANSLT
!
! Append appropriate coefficients (including decentering epsilon values -- See
! Notes 13,14,15)
!
  dt   = 0.5_r8*ztodt
  tmp1 = 0.5_r8*ztodt/coslat
  do k = 1,plev
     tmp = cappa*t0(k)*hypi(plevp)/hypm(k)
     do i = 1,nlon
        urhsl (i,k) = tmp1*urhsl (i,k)
        vrhsl (i,k) = tmp1*vrhsl (i,k)
        trhsl (i,k) = dt  *trhsl (i,k)
        tsld0a(i,k) = dt  *tsld0a(i,k)
        prhsl (i,k) = dt  *prhsl (i,k)
        t2    (i,k) = ztodt*t2   (i,k)
        fu    (i,k) = ztodt*fu   (i,k)
        fv    (i,k) = ztodt*fv   (i,k)
!
! Combine terms.
! (Time split the midpoint results between arrival and departure
!  points)
!
        tmp2          = lpsstar(i)*hortalc(k)
        u3m1    (i,k) = (u3  (i,k)*coslat + onemeps*usldm(i,k)*dt)/coslat + &
                                            onemeps*urhsl(i,k) + &
                                            onemeps*fu   (i,k)*0.5_r8
        v3m1    (i,k) = (v3  (i,k)*coslat + onemeps*vsldm(i,k)*dt)/coslat + &
                                            onemeps*vrhsl(i,k) + &
                                            onemeps*fv   (i,k)*0.5_r8
        t3m1    (i,k) = t3  (i,k) +         onemeps*tsldm(i,k)*dt + &
                                            onemeps*trhsl(i,k) - tmp2 + &
                                            onemeps*t2(i,k)*0.5_r8
        prhssld (i,k) = (psldm(i,k)*dt + prhsl(i,k))*onemeps
        tarrsld (i,k) = (trhsl(i,k) + &
                         hybm(k)*tmp*tsld0a(i,k))*onepeps + tmp2 + onepeps*t2(i,k)*0.5_r8
        parrsld (i,k) = (prhsl(i,k) - hybd(k)*tsld0a(i,k))*onepeps
        lnpssld (i,k) = logps(i) - lpsstar(i) - tsld0a(i,k)*onemeps
        fu      (i,k) = urhsl(i,k)*coslat*onepeps + fu(i,k)*coslat*onepeps*0.5_r8
        fv      (i,k) = vrhsl(i,k)*coslat*onepeps + fv(i,k)*coslat*onepeps*0.5_r8
     end do
  end do
!
  return
end subroutine grmult_run

!
!-----------------------------------------------------------------------
!

end module grmult
